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Abstract 

Fluctuations in the local Newtonian gravitational field present a limit to high precision mea- 
surements, including searches for gravitational waves using laser interferometers. In this work, we 
present a model of this perturbing gravitational field and evaluate schemes to mitigate the effect 
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by estimating and subtracting it from the interferometer data stream. The resulting improvement 
in the detector sensitivity should allow, for the first time, terrestrial interferometers to detect 
£N| gravitational radiation from cosmological sources. 

PACS numbers: 04.80.Nn, 95.55.Ym, 07.60.Ly, 42.62.Eh, 04.80.-y 
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I. INTRODUCTION 



Gravitational waves (GW) from astrophysical sources have the promise of revealing a 
rich new vision of the universe [1]. In the past decade, several kilometer sized terrestrial 
detectors of gravitational waves (such as TAMA300 [2], GEO600 [3], Virgo [4], and the Laser 
Interferometer Gravitational- wave Obervatory (LIGO) [5]), have come online and made 
searches for signals in the 50-10000 Hz band. The reach of these ground based detectors at 
low frequencies is limited by seismic and gravitational perturbations on the earth. A set of 
space missions (NGO [6], DECIGO [7]) are being pursued to search for signals in the 10~ 4 -1 
Hz band. 

The multi-stage vibration isolation systems [8, 9] developed for GW detectors should, in 
principle, be capable of reducing the direct influence of the ambient seismic noise to below 
the quantum and thermodynamic limits of the interferometers. Unfortunately, there is no 
known way to shield the detectors' test masses from fluctuating gravitational forces. As 
shown in Figure 1, our calculations estimate that the fluctuations in the local Newtonian 
gravitational field will be the dominant source of the mirror's positional fluctuations near 
10 Hz. This noise source has been referred to as gravity gradient noise or Newtonian noise 
(NN) in previous literature. 

Early estimates of NN by Weiss [10], Saulson [11], and Hughes and Thorne [12] have made 
increasingly better estimates of the seismic environment and thereby, the gravitational noise. 
In this work, we update the estimations of Newtonian noise as well as describing a means 
to subtract its influence from the data stream. 

II. SIMULATION OF SEISMIC NEWTONIAN NOISE 

Since Newtonian noise cannot, at this time, be directly measured, we must base our 
estimates of subtraction capabilities on simulated noise. We attempt to obtain a sufficiently 
accurate estimate of the NN based on information about its source, which, in this case, 
is the seismic field. Other contributions to NN such as building vibrations or air pressure 
fluctuations are not considered here because, as discussed in Appendix A, they are not 
expected to limit GW interferometer sensitivities in the next ~10 years. 

In this Section, we give a description of the time-series generator for the seismic fields, the 
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FIG. 1. (Color online). Strain noise spectral density of Advanced LIGO (BLACK). The sensitivity 
of the initial LIGO (PINK - dashed) is shown for comparison. The Newtonian gravity noise 
(GREEN) is dominating the Advanced LIGO noise near 10 Hz. Other traces shown are other, 
non-gravitational, limits to the sensitivity: direct seismic vibrations (BROWN), quantum radiation 
pressure and shot noise (PURPLE), mirror thermal noise (RED), and mirror suspension thermal 
noise (BLUE). 

associated NN, and other instrumental noise of the interferometer and seismometers. We 
also discuss the suitability of our simulation as an estimate of NN at the LIGO sites. The 
problem is set up as a full time-domain simulation of seismic fields and instrumental noise. 
Instrumental noise such as seismometer noise or test-mass displacement not generated by 
NN is treated as stationary. In contrast, we do not assume stationarity of the seismic field. 
Here the attempt is to simulate a seismic field that is comparable to previously measured 
seismic noise, and to make the content of the field as complex as possible in order to test 
NN subtraction schemes on challenging scenarios. Still, due to computational limitations, 
simplifications are necessary. In the general case, if ground displacement £(r,t) is weak, 
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then the test-mass acceleration due to NN can be estimated by the integral over the entire 
ground medium 

5a N N(ro, t) = G J dV r^rp (t(f, *) ~ 3(e r • £(f, (1) 

where po is the density of the ground, G is the gravitational constant, r the position of the 
test mass, r points to locations in the ground, and e r is the unit vector pointing from r to 
f [13]. This Equation is valid for arbitrary seismic fields and represents the noise imprinted 
on the test-mass due to NN. In our simulation, we only consider seismic fields composed of 
surface waves. This simplification is enforced by computational limitations since generating 
NN time series from simulated 3D seismic fields would require months-long simulation runs. 
We expect this assumption to be reasonable, since surface waves are expected to have much 
larger amplitudes than body waves [12], and so surface waves give the dominant contribution 
to NN at the surface, however seismic array measurements currently in progress at the LIGO 
sites will confirm this and several other of our assumptions for our particular sites. 

Freely propagating surface waves like Rayleigh waves and their overtones produce NN in 
such a way that there is always an effective 2D representation of the problem (which is not 
generally true for all supported wave fields, such as scattered waves). This implies that the 
numerical simulation can be set up as a 2D simulation, which is why Equations 2 and 3 only 
describe vertical displacement. This approach was chosen to reduce computational costs, 
and does not over simplify the subtraction problem as long as scattering of seismic waves is 
weak. 

The simulated seismic field is composed of two wave types; wavelets and symmetric 
surface waves. Wavelets represent seismic waves from far-field sources, while symmetric 
surface waves represent disturbances due to local sources. The vertical displacement due to 
a wavelet is described by 

£(f, t) = £o exp(-r 2 /(2AT) 2 ) cos(27r/r + fa) (2) 

with r = t — k ■ (r — r* )/(27r/). Twenty wavelets are injected for each second of time 
series. Frequencies / are drawn from a uniform logarithmic distribution between 8 Hz and 
30 Hz, which includes the full range of frequencies for which NN is expected to be dominant. 
Wavelet durations AT are uniformly distributed between 10/ / to 20//, to represent that 
wavelets can vary in duration depending on the type and source of the disturbance. The 
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distribution of wave vectors k is isotropic, to represent far-field sources from all directions. 
The average speed of sound for seismic waves in the ground is approximately 200 m/s [14] 
so we allow seismic speeds in our simulation to vary uniformly from 195 m/s to 205 m/s. 
This variance in speed is a brute-force method to simulate wave scattering, but it is very 
likely an overestimation of the effect. The initial location r*o of the wave maximum lies in 
the direction of the back-azimuth of the incident wave such that the wave is guaranteed to 
reach the location of the test mass within the simulated duration of the time series. The 
initial phase 0o of any single seismic wave is not a critical parameter. What is important 
is that not all seismic fields in our simulation have the same phase, so 0o is drawn from a 
uniform distribution between and 2ir. 

The second type of wave in the simulation is the symmetric surface wave, with vertical 
displacement 



with R — \r — Tq\. Equation 3 represents fields from sources located at r*o with distance tq 
drawn uniformly between 10 m and 20 m. All other parameters are obtained in the same way 
as for the wavelet, where as before the variation in seismic speed leads to a corresponding 
variation of the wave number fc . We assume that local sources do not vary strongly over the 
relevant time scales (defined by the subtraction procedure; see following Sections), so that 
the local sources are considered stationary. A fixed number of 10 waves from local sources 
is used. We believe that this combination of twenty wavelets and ten spherical symmetric 
waves is more complex than any typical seismic field we will see, but as with the relative 
amplitudes of surface vs. body waves, our array measurements will confirm our assumption. 

The simulation covers a surface area of 100 m x 100 m with the test mass at its center, 
which is larger than the area from which interesting NN contributions are expected [13]. 
The number of grid points along each dimension is N = 201 so that the grid-point spacing 
is 0.5 m. We choose a 201 x 201 point grid as a compromise between overall grid area, grid 
spacing, and computational time. The test-mass is suspended 1.5 m above ground, which 
is approximately the height of the LIGO test-masses. As the effective 2D representation is 
based on the surface term of the gravity perturbations and not the full dipole form [13], 
we convert the integral in Equation 1 into a discrete sum over grid nodes. Using only the 
surface contribution to the integral, the test-mass acceleration along the direction of the 
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interferometer arm is 
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(t) = G Po dS ]T 



Y~ cos (4>i) 



(4) 



i=l 



where dS is the area of the grid, £j(t) is the vertical displacement of grid point i at time t, 
is the distance between the grid point to the test mass, and fa is the angle between vector 
pointing from the test mass to the grid point and the direction of the interferometer arm. 
The sum over grid points in Equation 4 is used to determine the time series of the NN at the 
test-mass. Time series for each seismometer in Sections IV and V are calculated separately 
using Equations 2 and 3, so seismometer locations are not restricted to coinciding with grid 
points. 

We utilize models of the instrumental noise of seismometers and the strain noise of an 
interferometer to more accurately determine the NN subtraction efficacy, as described in 
Sections IV and V. The instrumental noise of all seismometers is simulated with spectral 
densities that are white (frequency independent) in units of velocity and have a value of 
10~ 10 m/v / Hz at 10 Hz. This is a conservative estimate for commercial geophones. 

The seismic spectrum itself plays a minor role for the purpose of this paper, but never- 
theless we defined distributions for £o in Equations 2 and 3 in such a way that the spectral 
density approximates the median spectrum measured at the LIGO sites. The plot in Fig- 
ure 2 shows the histogram of unaveraged 128 s spectra measured at the LIGO Livingston 
site over a time of one year during the last science run, and the black curve represents the 
average spectrum of the simulated seismic field. The average spectral density derived from 
the histogram is about a factor of 2 to 3 larger at frequencies between 10 Hz and 30 Hz than 
the model used in [12] with correspondingly larger NN spectrum. 

The strain noise model (excluding NN) that we use to simulate interferometer noise is: 



To convert from the displacement noise of a single mass into strain noise, we multiply by 2 
to account for the incoherent sum of 4 masses and then divide by the interferometer length, 
L = 4 km, to get strain. This is a representative noise model for proposed upgrades to 
Advanced LIGO and Advanced VIRGO, which we refer to as 3 rd generation ground-based 
detectors [15]. Future detectors built at new sites, such as the proposed Einstein Telescope, 





6 





10 5 

■ yj 












10 


> 








I 


.7 




10 


E 




3 




v_ 




eel 


10 8 

1 u 


SL 




W 




o 




ismi 


10" 9 











10 



10 




20 30 
Frequency [Hz] 



FIG. 2. (Color online). Histogram of one year of unaveraged 128 s seismic spectra measured during 
the last LIGO science run inside the corner station of the Livingston detector. The black curve is 
the spectral density of the simulated seismic field. The spectral histogram of the Hanford site is 
very similar for the frequencies plotted here. 



we call 4 th generation detectors [16]. 

The sampling frequency for all time series is / s = 100 Hz and the observation time is 
T = 100 s. We plan to test our subtraction techniques on longer duration simulated data in 
the future, however computational time restraints have kept us to this moderate duration 
for the time being. All time series are high-passed with corner frequency 5 Hz directly after 
being generated to avoid numerical problems. As can be seen in Figure 3, interferometer 
noise dominates NN below 8 Hz and above ~20 Hz so that we can safely ignore frequencies 
outside this range when testing NN subtraction methods. 
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FIG. 3. (Color online). Spectrum of simulated Newtonian noise (RED), reference 3 rd generation 
sensitivity curve (BLUE), and simulated interferometer noise based on this noise model (GREEN). 

III. SENSOR ARRAY OPTIMIZATION 



Optimization of seismic arrays with respect to NN subtraction was discussed in [17]. The 
authors calculated subtraction residuals analytically by evaluating explicitly the correlation 
between seismometers and the test-mass as a function of seismometer locations. The average 
subtraction residual can be written as 



R 



j _ ^SN ' (Css) 1 • ^SN 



C 



NN 



(6) 



Here, C*sn is the cross-correlation vector between seismometers and the NN acceleration of 
the test-mass, Css is the cross-correlation matrix between seismometers, and Cnn is the NN 
variance. These quantities can also be interpreted as (cross-) correlation spectral densities. 
Given a fixed number of seismometers, the optimal array is found by changing seismometer 
locations and minimizing \/R. The equation is idealized as it does not depend on any 
details about the way subtraction is implemented, i. e. whether a finite impulse response 
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FIG. 4. Comparison between the theoretical model for seismic correlation of isotropic plane- wave 
surface fields as described by Equation 7 (solid line), and the correlation calculated by the simu- 
lation of a field composed of wavelets and locally generated waves (dotted line). The correlations 
match at close distances and deviate strongly at larger distances. Since seismic fields in the con- 
text of NN subtraction only matter very near the test mass, the match between simulated and 
theoretical correlations at small distances means that the optimal array determined analytically by 
minimizing y/R in Equation 6 should also perform well in simulation, and more importantly that 
the simulation should be representative of our real subtraction ability. 

(FIR) filter is used, or some non-causal post-subtraction filter (see following two Sections 
for details). For this reason it describes the performance of all subtraction methods that are 
based on linear filtering, and the optimal array found by minimizing \/R is universal for all 
linear noise filters. Since it is very likely that different noise cancellation techniques will be 
combined in practice, it seems that optimization based on Equation 6 is the best one can 
do. 

Correlation patterns of surface waves observed in nature are often well approximated by 



Bessel functions that characterize isotropic plane- wave surface fields [18, 19]. Adopting a 
more convenient normalization the corresponding seismic correlation Css between two points 
r*j, fj on the surface is given by 

C ss (ri, fj) = Jo(27r|fi - f^/X) + ^2 (7) 

where A is the length of the seismic wave, and SNR is the signal-to-noise ratio of the seis- 
mometers. To find out how well this theoretical model approximates the seismic correlation 
in the simulation, we calculated Css between seismometers of increasing distance using our 
simulated seismic fields. The result is shown in Figure 4, where we show that the correlation 
vs. distance of our simulated seismic fields match the theoretical correlation of seismic fields 
fairly well, albeit not precisely. 

Other terms in Equation 6, using the same normalization as for Equation 7, are the NN 
variance 

C NN = 0.5 (8) 

and the correlation between seismic displacement and NN acceleration of the test-mass 
located at the origin 

Cm{n) = UZirn/Xp (9) 

with Ti = \fi\, and Xi is the projection of r*j onto the direction of the interferometer arm. 
Since Equation 6 is independent of seismic or NN amplitudes, we can use any suitable 
normalization of the seismic field or the NN. 

Finding the optimal array is not a trivial task. The result of a step-wise optimization by 
placing one seismometer after another leads to array configurations very different from the 
optimum. For the model described by Equations 7 to 9, the step-wise optimization yields 
a straight line of seismometers along the direction of the arm, approximately symmetric 
about the test-mass, independent of the number of seismometers. Therefore, configurations 
close to the optimum can only be found by optimizing all seismometer locations simultane- 
ously. A systematic numerical search for the optimum for more than a few seismometers is 
prohibitively computationally expensive, and approximate numerical optimization methods 
need to be applied. The array configuration that we call optimal in the following sections 
is shown in Figure 5. It was found numerically by running a particle-swarm minimization 
code [20, 21] to optimize the location of 10 noiseless seismometers. It should be clear that 
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FIG. 5. Locations of 10 sensors resulting from numerically minimizing the subtraction residual. 
The optimal array should be symmetric, but the subtraction residual is less than 10 -6 at 10 Hz for 
this array. The colors indicate the normalized seismic correlation between seismometer 1 and all 
other seismometers. 



the optimal array should have some kind of symmetry, so we know that this configuration 
is sub-optimal. The optimization was stopped at a residual y/R ~ 10~ 6 at 10 Hz, which is 
negligible for all practical purposes. So while this configuration does not represent a global 
optimum, its subtraction performance should be sufficient for Advanced LIGO and 3 rd gen- 
eration detectors. As many configurations yield similarly small subtraction residuals, we 
added further components to the cost function \/R to make sure that seismometers are not 
placed too close to each other. The array shown in Figure 5 is the result of minimizing this 
combined cost function. 

As one can see from Equations 7 to 9, the residual R is a function of seismic wavelength, 
and therefore frequency, and broadband subtraction performance needs to be investigated. 
The subtraction residual of the array in Figure 5 was minimized at 10 Hz for a seismic 
wave speed of 200 m/s. In Figure 6 we show the subtraction residual as a function of 
frequency for various array configurations. One can see how the number of seismometers 
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FIG. 6. Subtraction residual as denned in Equation 6 vs. frequency for the array shown in 
Figure 5, and three different spiral configurations. It is assumed that the seismometers measure 
ground motion with SNR = 100 at all frequencies. The Rayleigh-wave speed is 200 m/s. 



and the array size affect subtraction residuals. It is clear that a very small array does 
not perform well at low frequencies since it provides highly degenerate information at these 
frequencies whereas larger arrays sample a larger part of the seismic wave. A smaller number 
of seismometers simply leads to a broadband increase of subtraction residuals except for the 
smallest frequencies. We want to emphasize that these theoretical predictions only hold 
approximately for the numerical simulation presented in the following sections, since it does 
not account for details of the subtraction method as explained before. 



IV. OFF-LINE POST SUBTRACTION 

For the purpose of this paper, off-line post-subtraction denotes the cancellation of noise 
in recorded data. The noise cancellation filter can therefore be causal or non-causal. In 
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this section, we will present a simple non-causal implementation of the post-subtraction. 
The method chosen here is to cancel NN on short segments of recorded interferometer 
time series. The basic idea is to optimally construct a vector of filter coefficients, one 
coefficient per seismometer, and then use these coefficients to form a linear superposition of 
the seismometer channels as the NN estimate. For a general introduction to digital filtering 
techniques, please see, for example, [22]. 



The residual R in Equation 10 corresponds to the interferometer data / minus the NN 
estimate from the seismometer data. The time series used to calculate the NN estimate are 
pre-conditioned with whitening and band-pass filters focussing on the 8 Hz to 30 Hz NN band 
before the correlations are evaluated. We also found it necessary to apply an anti-aliasing 
window (we used the high-gain Nuttall window) for reasons that will be described below. 
All quantities subject to the preconditioning are marked with a "c". (I c ,S c )i denotes the 
vector of cross correlations between the interferometer data P and all seismometers S c using 
data of segment i acquired between tj and t i+ i. Similarly, (S c ,S c )i is the cross-correlation 
matrix between all seismometers. This means that the filter used here will have one filter 
coefficient per seismometer for the entire time interval ti to U+i. 

We must determine a reasonable time duration for each segment. Segments are too short 
if the spectral resolution is too small to disentangle seismic waves at different frequencies. 
Segments may be too long if the number of seismic waves in that time frame becomes 
large. A Wiener filter that sees many seismic waves may begin to average over the different 
waveforms and provide non-optimal noise suppression. Choosing the goldilocks segment 
duration is somewhat arbitrary, however it is likely that the appropriate duration depends 
as much on the nature of the seismic field as on the frequency band targeted by the filter. 
With our simulation we found the best subtraction performance for 2 s long segments. This 
is an acausal technique, so testing can be done offline to determine the duration for which 
we see maximal NN suppression on the real data. 

Since filter coefficients are re-evaluated for each segment i, a simple subtraction of NN 
estimates from consecutive segments can lead to discontinuities in the residual time series. 
For this reason the Nuttall anti-aliasing window is applied so that noise subtraction is 
suppressed at the beginning and end of a time segment. Consequently time segments are 
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defined with overlap to provide continuous subtraction of NN. Using the Nuttall window, we 
found excellent subtraction performance with 0.3 fractional segment overlap. Again, some 
investigation can be done to optimize this number for real data in the future. 

Optimal array design has already been discussed in Section III. We will compare the 
subtraction performance of the optimal array presented there with a circular, a spiral, and 
a linear array. All arrays contain 10 seismometers, and are optimized in terms of the extent 
of the array relative to the location of the test-mass. The linear array is simply a line 
of uniformly spaced seismometers along the direction of the arm extending 8 m away from 
the test mass in both directions. This linear array is slightly different from the result of the 
stepwise optimization discussed in Section III, but the subtraction residuals are similar. The 
circular array consists of one seismometer under the test mass and 9 seismometers in a circle 
of radius 5 m around the test mass. The configuration of the spiral array is shown in Figure 8 
of the following section. The residuals of the noise subtraction (described in Equation 10) 
for each array are shown in Figure 7. The noise model represents the sensitivity curve 
of a potential upgrade of the advanced detectors not including the NN, as described by 
Equation 5. Approximately, all arrays perform equally well in post-subtraction. The goal to 
reduce the NN residuals to a level below the noise model is achieved over the entire NN band 
except for the very smallest frequencies. In Section V, we will investigate the possibility of 
combining the post-subtraction with an online feed-forward cancellation. 

V. ONLINE FEED-FORWARD SUBTRACTION 

Online feed-forward (FF) subtraction can be implemented in two ways. It is possible to 
continuously cancel NN by exerting a cancellation force directly on the test-masses. Alter- 
natively, the cancellation can be done on interferometer data. If we had ideal, noise-free 
actuators, the residuals resulting from applying forces on the test-masses and online FF 
cancellation applied to the data would be the same. Applying hardware cancellation forces 
could also be used to suppress the problem of any non-linear response of the detector to 
strong NN forces, but is very technically challenging to implement [23]. Since we do not 
believe that near-future detectors will suffer from significant non-linear upconversion due to 
NN, we only consider FF cancellation applied to the data. 

The main difference between online FF and post-subtraction is that online subtraction 
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FIG. 7. (Color online). Offline Newtonian noise subtraction efficacy for an upgrade of the advanced 
detectors. Spectrum of simulated Newtonian noise (GREEN), proposed 3 rd generation sensitiv- 
ity curve (BLUE), and NN residuals of post-subtraction for a spiral array (RED), circular array 
(CYAN), linear array (MAGENTA), and the optimal array (BEIGE). Filters derived from all four 
arrays reduce the simulated NN to a level below other sources of interferometer noise as represented 
by the noise model. 




can only be done with causal niters. Furthermore, the FF filter coefficients can only change 
slowly in time following slow changes of average correlations between seismometers and the 
test mass. 

The FF subtraction scheme that we propose is based on a multi-input, single-output 
(MISO) finite-impulse response (FIR) filter that is continuously applied to the interferom- 
eter output to filter out the NN as was already demonstrated successfully for seismic noise 
cancellation schemes [24]. The inputs consist of the seismometer channels, and the single 
filter output is the NN estimate. 

Average correlation between seismometers and interferometer data has a predictable form 
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since average properties of the seismic field depend solely on the wave composition of the 
seismic field, which is characteristic for each site. This correlation pattern was investigated 
in Section III, where we showed that the simplest theoretical model is a good representation 
even for the more complex wave composition that is used in our numerical simulation. 

As we will show in the following, sufficient FF subtraction down to the level of other noise 
contributions can be achieved with a variety of array configurations including arrays that 
have seismometers with negligible correlation with the test mass NN. The more important 
design factors are the number of seismometers and the size of the area covered by the array. 

The only filter parameter that is predefined is the order of the FIR filter, i. e. the number 
of filter coefficients. The filter order essentially determines the time span of the filter. 
Therefore, similar to the post-subtraction scheme, we found that the order can be too high, 
in which case the seismic array cannot provide sufficient information to disentangle NN 
contributions from individual seismic waves. The filter order is too low when an insufficient 
amount of data is used to accurately estimate the NN from individual, resolved waves. We 
will later explain why the wave nature of the seismic field still matters in the context of 
FF cancellation. The FIR filter that yielded sufficient subtraction in all simulation runs has 
order N = 50 corresponding to a time span of 0.5 s. The MISO FIR filter coefficients were 
calculated from the 100 s long seismometer and test-mass time series generated as described 
in Section II. All time series are pre-conditioned with band-pass and whitening filters. An 
example of a Bode plot of the filter for a spiral array is shown in Figure 8. The fact that for 
example seismometers 3 or 5 have relatively high filter magnitudes at some frequencies is 
interesting since their correlation with the NN is very small (as calculated by Equation 9). 
This situation can be described as a trade-off between gaining information about how NN 
is generated close to each seismometer (the simple local model), and gaining information 
about how NN integrates over the seismic field based on its wave nature. 

The FF noise cancellation performance is shown in Figure 9. Since the FIR filter coeffi- 
cients are the same for the entire time series (see Section VI for alternative filter implemen- 
tations), we included two NN residuals, one for the Wiener filter that subtracts on the same 
time series used to calculate the filter coefficients, and a second one where the same filter 
is applied to subtract NN from another time series. The two time series represent different 
sets of local sources and wavelets. The subtraction performance is very similar for the two 
cases, and therefore we can conclude that subtraction performance does not depend as much 
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on the specific wave content of the seismic field as it depends on the average correlations 
between sensors and the NN. While the Wiener filter applied to the data on which it was 
trained use of the filter and could not be applied online, it is useful to see that 

the subtraction efficacy does not degrade for times that are not the training data for the 
filter. As with the post-subtraction, FF cancellation performed similarly for the circular, 
linear and spiral arrays. 

Finally, we investigate the possibility of combining the online FF cancellation with post- 
subtraction. Figure 10 shows the residual NN spectra for the three subtraction methods. 
Overall, there is no clear advantage or disadvantage to combining the two methods. When 
both techniques are applied, NN residuals are smaller at lower frequencies, but residuals 
are larger at higher frequencies. In conclusion, it was demonstrated that the standard 
static MISO FIR Wiener filter provides robust and sufficient subtraction results. Whereas 
a combination of FF and post-subtraction does not give further improvement in simulation, 
it could prove more effective in scenarios where strong occasional seismic disturbances leave 
significant residuals after FF cancellation. 

VI. FUTURE SIMULATION WORK 

In this section, we outline the main challenges associated with NN subtraction and discuss 
in more detail where our numerical simulation needs to be refined to better understand the 
associated risks. By far the greatest challenge of NN subtraction is to make sure that all 
relevant sources of gravity perturbations are identified. In this context, our simulation is 
certainly highly simplified. However, our estimates indicate that the seismic field gives the 
only NN contribution that will be relevant to the advanced detectors or their potential 
upgrades as represented by the strain-noise model used in this paper (see Appendix A for 
details). This justifies the exclusion of other NN sources in our numerical simulation. 

The more interesting aspect of the source-identification problem is whether we can be 
sure that all relevant degrees of freedom of the seismic field are or can be monitored. The 
seismic array in Section III is designed based on prior knowledge. For example, in this 
simulation we assumed that there are no significant NN contributions from body waves that 
propagate through the ground in all directions. Surface waves are expected to have much 
larger amplitudes than body waves near 10 Hz at the LIGO sites [12], but a detailed study 
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of the fields has not yet been done. Our measured seismic spectra and NN estimates shown 
in Figure if do not tell us the wave content of the seismic field. In the near future, we 
will construct a seismic array at the LIGO Hanford site from which we should be able to 
determine the relative significance of body waves. 

Another related issue that is often discussed is the scattering of seismic waves, which 
we assume is negligible for our simulation. This could indeed pose a severe problem to 
NN subtraction even if scattering is identified and fully characterized. Scattering can in 
principle make it impossible to estimate NN from seismic measurements at the surface 
since it can lead to a more complex field structure that is not completely characterized by 
surface displacement. Moreover, it is possible that scattered waves have higher wave numbers 
compared to the freely propagating surface waves, so that the density of the seismic array 
would need to be increased to a point where it becomes infeasible or at least very challenging 
to monitor the entire field accurately. However, for the modest subtraction performance 
that we require for the detector model we consider in this paper, the major portion of the 
gravitational noise perturbation comes from the surface area very close to the test-mass, as 
opposed to higher subtraction factors, where a substantially larger surface area needs to be 
considered. It follows that scattering will only be a problem if it is strong enough to alter 
seismic waveforms significantly over very short propagation distances. Since the ground 
medium close to the test mass is fairly uniform, high scattering cross sections are unlikely 
to be observed. Seismic measurements will be necessary to test this assumption. 

Methods that we have not investigated in this paper could help to mitigate some of the 
risks if necessary. In our simulation the FF filter used was implemented as a static Wiener 
filter, however it is possible to let the filter coefficients adapt slowly to changes of the seismic 
field. This adaptive filter technology has many applications and is well established [25]. 
Also, once the array design has been chosen based on previous seismic measurements, cross- 
correlations observed with this array can help to find better array configurations. In other 
words, it will be possible to adapt to changing properties of the seismic field not only through 
adaptive filter technologies, but also through changes in the hardware configuration. 

Adding more details to the numerical simulation like scattered waves or body waves 
can tell us in advance how the array would have to be modified to maintain the same 
level of subtraction performance for these more complicated scenarios. It would be wise to 
investigate all scenarios irrespective of how well we think we understand the seismic fields. 

18 



VII. CONCLUSIONS 



We have shown that a relatively small number of medium sensitivity geophones or ac- 
celerometers can be used to estimate the Newtonian gravitational fluctuations with a rea- 
sonably high accuracy. Under our simplifying assumptions for the seismic fields and the 
structure of the ground, this allows us to use seismic data to subtract the gravitational noise 
due to seismic motion from the interferometer data stream well enough that the Advanced 
LIGO, as well as 3 rd generation detectors should not be limited by this terrestrial noise 
source. 

We found that the array configuration has a minor impact on the subtraction residuals. 
The more important design parameters are the number of seismometers, the area covered 
by the seismic array, and proper pre-conditioning of the time series that are used for the NN 
estimate. 

Our numerical simulation needs to be developed further to test subtraction of other 
possible contributions to the seismic field that have mostly been considered insignificant 
for the NN problem in advanced detectors in the past, as for example body waves and 
scattered waves. Testing cancellation of NN by factors of 10 or more requires a more accurate 
simulation of seismic fields. 

The offline, acausal subtraction scheme should naturally outperform the online, adaptive 
causal feed-forward technique, but for the simple implementation of the post-subtraction 
used in this paper, the subtraction performances were comparable. To get latency for a 
cleaned-up data stream to be less than ~1 minute, we will do initial subtraction online and 
then make the final subtraction offline. 

These NN subtraction techniques will have a modest improvement on the the second 
generation detectors (Advanced LIGO, Advanced Virgo, KAGRA), but the true promise 
will come towards the end of the decade. At that time these techniques will be necessary 
to achieve the next order of magnitude improvement in astrophysical reach with third- 
generation of detectors. 
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Appendix A: Newtonian-noise budget for the LIGO sites 

In 2011, several measurements were carried out at the LIGO sites to define a Newtonian 
noise budget [26]. Using accelerometers, vibrations were monitored on water pipes, near 
exhaust fans, on top of the buildings and on the walls. With microphones, sound spectra 
were measured inside and outside of the LIGO buildings. The resulting NN estimates for 
each of these sources are summarized in Figure 11. Note that the seismic curves for both 
LIGO sites presented in Figure 11 are more recent, using a more accurate, non-averaged, 
analysis of the seismic percentiles, as compared to [26]. In addition, the plot contains the 
3rd generation strain noise model and estimates of seismic NN for both sites. All curves are 
derived from 90th percentiles of spectral histograms similar to the one shown in Figure 2 for 
the seismic measurement at the LIGO Livingston site. According to these estimates, seismic 
NN is the only significant noise contribution for 3 rd generation and earlier detectors. 
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FIG. 8. (Color online). The upper plot shows the configuration of the spiral array. The colors 
correspond to the normalized seismic correlation between all seismometers and seismometer 1. The 
numbering of seismometers corresponds to the traces in the lower plot, which shows the magnitude 
of the FIR filter for each sensor in units of test-26ass (NN) displacement over seismic displacement. 
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FIG. 9. (Color online). Spectrum of simulated Newtonian noise (GREEN), proposed 3 rd generation 
sensitivity curve (BLUE), and NN residuals of FF subtraction on the training set (RED), and on 
a second set of time series using the same filter (CYAN) using the 10 sensor optimal array. 
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FIG. 10. (Color online). Feed- forward Newtonian noise subtraction efficacy for 3 rd generation 
LIGO detectors. Simulated NN before subtraction (GREEN), expected strain sensitivity (BLUE), 
NN residuals after subtraction using post-processing (RED), online feed forward (CYAN), and both 
methods combined (MAGENTA). Note that the combination of methods is close to the same level 
as either method individually. This indicates that we can safely apply feed forward subtraction in 
realtime, and clean up any leftover noise in post-processing if needed. 
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FIG. 11. (Color online). Seismic NN estimates for the LIGO Livingston (LLO) and Hanford 
(LHO) sites (GREEN and RED), the 3 rd generation strain noise model (BLUE), and additional NN 
estimates from vibrations of walls (CYAN), building tilts (MAGENTA), exhaust fans (BEIGE), and 
sound waves inside buildings (BLACK). Seismic contributions are the only NN source significant 
for 3 rd generation detectors and earlier. Building tilt will be important for detectors beyond the 
3 rd generation, but is not a dominating noise source at this time. 
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